Streszczenie raportu

Raport powstał na bazie zbioru ponad 50 tysięcy rekordów z kilkudziesięciu ostatnich lat, dotyczących rozmiarów łowionych w tym okresie śledzi. Dane zostały uzupełnione o brakujące wartości, poddane analizie pod kątem korelacji poszczególnych zmiennych i wykorzystane do zbudowania przez uczenie maszynowe metodą Random Forest regresora pozwalającego przewidywać rozmiar złowionej ryby na podstawie parametrów środowiskowych.

Przeprowadzone badanie wykazało, że istnieje w szczególności silne powiązanie pomiędzy stopniowym wzrostem mierzonych temperatur przy powierzchni wody, a trendem spadkowym w długości łowionych śledzi.

Zastosowane biblioteki, czynności wstępne

Wczytanie bibliotek dplyr (przetwarzanie danych), ggplot2 i plotly (prezentacja graficzna), easyGgplot2 (rozszerzenia do ggplot2) oraz caret. Dodatkowo wykorzystanie doMC do przetwarzania równoległego.

library(dplyr)
library(plotly) # devtools::install_github(“ropensci/plotly”)
library(ggplot2)
library(easyGgplot2) # devtools::install_github("kassambara/easyGgplot2")
library(reshape2)
library(caret)
library(doMC)
registerDoMC(cores = 4)

Wczytanie i wstępnie przetworzenie danych

raw.data <- read.csv('sledzie.csv', na.strings = '?')

Dane zawierają sporadycznie występujące puste wartości, które zostaną wypełnione medianami występującymi dla odpowiednich wartości kolumn xmonth i recr, które nie zawierają wartości pustych. Uznano to rozwiązanie za akceptowalne dla celów tej pracy, ponieważ liczności zbiorów wartości poszczególnych kolumn są niewielkie (w granicach 50 unikalnych wartości) w porównaniu z liczbą wierszów całej tabeli (ponad 50 tysięcy). Brakujące wartości są na tyle nieliczne, że wpływ ewentualnych przekłamań na wynik pracy można traktować jako pomijalny.

processed.data <-
  raw.data %>%
    rowwise() %>%
    mutate_all(
      funs(
        replace(., is.na(.), (
            function(rownum, v, ref_recr, ref_xmonth) {
              if (!is.na(v)) {
                return(v);
              }
              else {
                filtered = raw.data %>% filter(recr == ref_recr, xmonth == ref_xmonth)
                med = (filtered %>% summarize(median(., na.rm = TRUE)))[[1]]
                return(med)
              }
            }
          ) (X, ., recr, xmonth)
        )
      )
    )

Podstawowe podsumowanie danych

Zbiór zawiera 52582 rekordów. Poniższa tabelka zawiera podstawowe statystyki dotyczące zbioru danych uzupełnionego o brakujące wartości komórek.

knitr::kable(summary(processed.data[, 1:8])); knitr::kable(summary(processed.data[, 9:16]))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2
Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808
Median :26290 Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750 Median :21.435 Median : 7.0000 Median :24.859
Mean :26290 Mean :25.3 Mean : 0.4440 Mean : 2.0257 Mean : 9.994 Mean :21.219 Mean : 12.7967 Mean :28.422
3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232
Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736
fbar recr cumf totaln sst sal xmonth nao
Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
Median :0.3320 Median : 421391 Median :0.23191 Median : 539558 Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
Mean :0.3304 Mean : 520366 Mean :0.22981 Mean : 514973 Mean :13.87 Mean :35.51 Mean : 7.258 Mean :-0.09236
3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000

Rozkłady wartości zmiennych

Wartości średnie oznaczone liniami przerywanymi w kolorze czerwonym.

plot.var <- function(varname, ...) {
  ggplot(processed.data, aes(x=get(varname))) +
    geom_histogram(color="black", fill="white", ...) +
    geom_vline(aes(xintercept=mean(get(varname))), color="red", linetype="dashed", size=1) +
    xlab(varname)
}

ggplot2.multiplot(
  plot.var("length", binwidth = 0.5), plot.var("cfin1", bins = 30), plot.var("cfin2", bins = 15),
  plot.var("chel1", bins = 30), plot.var("chel2", bins = 30), plot.var("lcop1", bins = 30),
  plot.var("lcop2", bins = 30), plot.var("fbar", binwidth = 0.08), plot.var("recr", bins = 15),
  plot.var("cumf", bins = 9), plot.var("totaln", bins = 12), plot.var("sst", binwidth = 0.15),
  plot.var("sal", binwidth = 0.012), plot.var("nao", binwidth = 1),# plot.var("xmonth", binwidth = 1),
  cols = 3
)

Korelacja zmiennych

cor_matrix = round(cor(processed.data), 2)
#cor_matrix[lower.tri(cor_matrix)] <- NA
melted_cor_matrix <- melt(cor_matrix, na.rm = TRUE)

ggplot(data = melted_cor_matrix, aes(Var2, Var1, fill = value)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "#91bfdb", high = "#fc8d59", mid = "#ffffbf", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Korelacja\nPearsona") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
        size = 12, hjust = 1)) +
  coord_fixed()

Mapa korelacji pozwala na wyciągnięcie pewnych wstępnych wniosków dotyczących przydatności poszczególnych zmiennych do budowy regresora przewidującego długość łowionych ryb (atrybut length).

processed.data <- processed.data %>% select(-lcop1, -lcop2, -cumf, -xmonth)

Ilustracja zmienności długości łowionych ryb

Poniższy wykres przedstawia zależność długości ryby (mierzonej z dokładnością ) od czasu (numeru sekwencyjnego próbki w zbiorze). Dodatkowo, próbki zostały pokolorowane w zależności od sst (temperatury przy powierzchni wody) - jaśniejszy odcień koloru niebieskiego oznacza wyższą temperaturę, a u dołu wykresu znajduje sie pasek pozwalający na przełączanie pomiędzy próbkami w momentach, gdy oscylacja północnoatlantycka (atrybut nao) osiągał wartości dodatnie, a tymi wykonanymi przy ujemnej wartości tej zmiennej (nao_range odpowiednio 1 i -1).

Na wykresie ze względu na wydajność nakreslono co 35 rekord ze zbioru.

ggplotly(
  ggplot(processed.data[seq(1, nrow(processed.data), freq), ] %>% mutate(nao_range = sign(nao)), aes(X, length)) +
    geom_point(aes(color=sst, frame=nao_range)) +
    labs(x = "numer sekwencyjny próbki",
         y = "długość ryby (cm)",
         color = "Temp. pow.\nwody (st. C)",
         title = "Długość ryby, czas i temperatura")
)

Wizualizacja zdaje się być kolejnym czynnikiem mogącym potwierdzać przypuszczenie o związku temperatury przy powierzchni wody z długością łowionych ryb. W początkowym okresie badania (około 17000 próbek) wartość temperatury wykazywała tendencję spadkową (zilustrowane ciemniejącym odcieniem błękitu punktów), rosły zaś wskazania pomiaru długości ryb. W dalszej części wykresu temperatura przy powierzchni wody rośnie, spadają zaś pomiary długości śledzi.

Trudniejszy do oszacowania wydaje się być wpływ parametru oscylacji północnoatlantyckiej. Jest to zjawisko cykliczne, która to cykliczność polega na okresowej dominacji odczytów dodatnich lub ujemnych - co zresztą znajduje potwierdzenie w wykresie, jeśli spojrzeć na okresowe zagęszczenie próbek dla ujemnego i dodatniego nao. Atrybut ten jest też w dość dużym stopniu skorelowany z sst - w okresie dominacji ujemnych odczytów nao dominują też niższe wartości sst i większe długości łowionych śledzi.

Konstrukcja regresora dla zmiennej length

Podział zbioru danych wejściowych

Podział zbioru w proporcjach: 60% zbiór uczący, 20% zbiór testowy, 20% zbiór walidacyjny.

set.seed(23)
inTraining <- createDataPartition(y = processed.data$length, p = 0.6, list = FALSE)
training <- processed.data[inTraining, ]
testing.and.validation <- processed.data[-inTraining, ]
inTesting <- createDataPartition(y = testing.and.validation$length, p = 0.5, list = FALSE)
testing = testing.and.validation[inTesting, ]
validation = testing.and.validation[-inTesting, ]

Schemat uczenia

Uczenie z wykorzystaniem wielokrotnej oceny krzyżowej.

ctrl <- trainControl(
  method = "repeatedcv",
  number = 2,
  repeats = 5
)

Uczenie metodą Random Forest

set.seed(23)
fit <- train(length ~ ., data = training, method = "rf", trControl = ctrl, verbose = TRUE, importance = TRUE)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
pred <- predict(fit, validation)
measures = postResample(pred = pred, obs = validation$length)
rmse <- measures[[1]]
rsquared <- measures[[2]]

Wygenerowany regresor dla zbioru walidującego osiągnął RMSE=1.1204929 i R2=0.5384306. Biorąc pod uwagę zakres wartości przewidywanej zmiennej length, taką wartość pierwiastka błędu średniokwadratowego można uznać za zadowalającą przynajmniej jako początkowy punkt odniesienia do analizy danych. W przypadku osiągniętej wartości R2 wyliczony model wyjaśnia około 50% wariancji zmiennej length i około 30% jej odchylenia standardowego. Nie jest to zatem predyktor zapewniający bardzo wysoką pewność otrzymania prawidłowego rezultatu, niemniej jednak warto wyciągnąć także wniosek z poniższego wykresu.

ggplot(validation %>% mutate(predicted_length = pred), aes(X, length)) +
  geom_point(aes(color = predicted_length - length)) +
  scale_colour_gradientn(colours = c("red","green",  "blue")) +
  labs(x = "numer sekwencyjny próbki",
     y = "długość ryby (cm)",
     color = "Błąd predykcji",
     title = "Relacja długości ryby do wartości przewidzianej przez model")

Wykres ten ponownie przedstawia zmienne w czasie zmierzone długości ryb, jednak tym razem punkty zostały pokolorowane stosownie do błędu predykcji dla danej próbki. Uzyskany model jest dosyć konserwatywny i niechętnie wychodzi poza pewien zakres wartości (na wykresie kolor zielony), notując w nim wysoką dokładność, zaś myląc się głównie w przypadkach skrajnych. Można przypuszczać, że liczność tych ostatnich jest główną przyczyną dosyć przeciętnego rezultatu współczynnika determinacji.

Ważność atrybutów

importance <- varImp(fit)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
plot(importance)

Analiza ważności poszczególnych zmiennych potwierdza wykorzystanie przez predyktor w największym stopniu zmiennej sst; druga co do bezwzględnej wartości korelacji zmienna nao nie została zaś w ogóle wykorzystana. Pozwala to na potwierdzenie przyjmowanego dotychczas przypuszczenia, że wzrost temperatury przy powierzchni wody ma kluczowe znaczenie dla obserwowanego trendu malejących wymiarów łowionych śledzi.